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ABSTRACT 


In an electrogasdynamic (EGD) generator the flow of a neutral medium 
carrying unipolar electrical charges sonstaeute: an electrical current. 
A model of the charge flow in an EGD generator was constructed for use 
in a computer program. The program is designed to solve Poisson's and 
Laplace's equations for both axisymmetric and two-dimensional geometries. 
Schlieren photographs of the charge cloud were used to determine the 
charge cloud profile required by the program. Computer generated pre- 
dictions agreed with three known solutions to Poisson's equation. 
Computer predictions of the effects of space charge flow modification 
were obtained. Space ahneae flow was modified both by increasing flow 
Speed and by manipulating the space charge electric field. Experimentally, 
this was accomplished by increasing flow stagnation pressure and by 
application of'‘a separate controllable electric field. Experimental 
results compare favorably with computer predictions. Some measurements 


were also made of the mobility range of the charged particles. 
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I. INTRODUCTION 


Electrostatic phenomena involve many everyday experiences, such as 
walking on a thick rug and then touching a grounded object. Practical 
applications of electrostatic devices range from Xerography to pollution 
control devices or to the van der Graaf generator [1]. Recently, interest 
in electrical power generation has brought attention to the electrogas- 
dynamic (EGD) generator. This direct energy conversion device is ana- 
logous to the van der Graaf generator and operates by convecting charges 
with a gas against an electric field to a point of collection. The 
charge flow constitutes an electrical current. The charge movement is 
accomplished by an exchange of momentum between the flowing dielectric 
medium and the charges. 

The common denominator of electrostatic devices is the existence of 
free charges in a neutral medium or dielectric -- for the present pur- 
poses a fluid. The flow of a medium with unipolar charges is termed 
Space charge flow and the understanding of this flow is critical to the 
design of the EGD generator. This thesis consists of a computer model 
of the space charge flow to be used as an EGD generator design tool. 

The computer model has some limitations inherent to the programming 
itself and these are discussed in Section III. In addition, parameters 
such as charge mobility and velocity as well as the inputs of geometry 
and jet configuration are needed. These are obtained with an experimental 
Set-up discussed in Section IV. The program is used to predict regions 
of breakdown, given an electrode geometry and a space charge cloud. In 
addition, the program is used to predict the effect of guard electrodes. 


These electrodes are placed for the purpose of ‘smoothing’ the electric 





field in the conversion region, and both computer predictions and 
experimental results are given in Section VI for the effect of the 


guard ring. 





II. THEORETICAL CONSIDERATIONS 


The EGD process requires producing, transporting and collecting the 
charges which comprise the electrical current. One means by which 
charges may be produced is by a corona discharge. Basically this scheme 
consists of two electrodes, which may be a needle and a concentric ring; 
the concentric ring or attractor electrode, when charged to a high posi- 
tive electrical potential with respect to the grounded needle, induces the 
corona discharge. The corona breaks down a neutral gas initially into 
its electrical components: molecular ions and free electrons. The corona 
needle will attract the positive ions, while negative charges will tend 
to migrate to the eerie ring. The EGD objective, instead, is to 
transport a current of negative charges away from this attractor ring. 

An ion's susceptibility to movement in an electric field is termed 
it's 'mobility', (k). Thus the drift velocity of the ion in the corona 
field is determined by the product of the ion mobility and the electric 
field strength to which the ion is exposed. Since the mobilities of 
molecular ions and free electrons are high, they would migrate rapidly 
to the corona electrodes and be lost to EGD power generation. However, 
when a supersaturated vapor is present, negatively charged droplet 
(micron) sized particles are nucleated which, for a given charge, will 
have fueh lower mobility than an ion [2,5]. Hence, the charged particle 
is more suited to EGD applications [6]. To move the charged particle 
against an electric field, the particle must be well coupled to the 
dielectric flow. When the dielectric gas transports the charge to a 
downstream collector, against an electric field (E), the total charge 


velocity (v,) may be described as the sum of the dielectric velocity 
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(v_) and the drift velocity (kE). 


fe + kE (1) 
In order to obtain the highest degree of coupling, the drift velocity 
of the charged particle must be low compared to the dielectric velocity 
{4]. The charged particles are carried out of the charge production 
section against an electric field in the conversion section, to the charge 
collector. This collector consists of an equipotential metallic surface. 
A path through a load toground is provided to complete the circuit 
started at the grounded corona electrode. See Figure (1) for a 
schematic of the EGD generator. 

There are several limitations that pertain to the conversion process. 
For an EGD generator to operate as intended, the neutral gas must remain 
non-conducting. This implies that the electric fields to which the 
dielectric is subjected must remain below the dielectric strength, or 
breakdown field strength, (E.)> of the gas. Fields above the dielectric 
strength will produce regions of high conductivity through ionization 
[3]. The dielectric gas or fluid is subjected to the field resulting 
from the applied potential through which the charged particles are carried 
and to internal fields generated by the presence of the charged particles. 
Breakdown from external fields occurs as a result of the electrode geo- 
metry. For a given voltage, the existence of a small radius of curva- 
ture will produce a large electric field. The field will exceed the 
dielectric strength when the radius of curvature is small enough. The 
range of allowable operating voltages is thus a function of the elctrode 
Spacing. In turn, the spacing of the electrodes is also a function of 
the fluid dynamics of the dielectric, of possible charge depletion 


mechanisms, of the operating pressures, and of generator size requirements. 
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Total charge per particle is limited by the balance between the 
particle surface tension and charge es enER forces, and 1s known as 
the Rayleigh limit. In practice, however, the total charge per particle 
is governed by the charge production method. For particles charged in 
a corona field, and which grow to micron size, the number of charges is 
approximately two orders of magnitude less than that predicted by the 
Rayleigh limit [6]. These limits must be considered in the modeling of 
the space charge flow in an EGD device. 

In general, the equations needed to describe the EGD flow will include 
the momentum, continuity and energy equations. These would be necessary 
to determine the dielectric jet profile and hence the space charge 
density distribution. However, if the profile can be determined 
experimentally, and if the energy drained from the flow is low, the 


required equations are reduced to three: 


Faraday's Law 


electric field (V/m) 


eee seo E = 
ot B = magnetic flux density vector (Wb /mé ) 
Gauss’ Law 
=— _ Pe 3 
V-E= = eo = Space charge density (C/m~) 
O Esha dielectric constant for vacuum and gas 
(V-m/A-s) 


Charge Conservation 
ORs Vp = 0 2s) 
where the total velocity of the charged particle is given by: 
Vy a iil K(E op + Ee) (3) 
and 1s a function of the dielectric flow speed, particle mobility, 


applied electric field eee and the space-charge-produced electric 


field (Eo). 


HZ 








For the EGD generator, B = 0, and Faraday's Law takes the form, 
yx E=0. This allows E to be defined as the negative gradient of some 
potential function 4, so that E = - Vd. Substituting into Gauss' Law, 


the resulting equation is known as Poisson's equation, 


a . (4) 


This equation is valid in the two regions of interest within the EGD 
channel. One, the neutral gas jet, contains the charged particles, while 
the second, outside the gas jet, contains no space charged. In the 
space-charge-free region, p, = 0, and Eq. (4) reduces to the more famil- 
iar Laplace equation. Figure (1) depicts the two regions. 

Although the kinetic energy of the injected gas jet is generally 
augmented by the transfer of momentum from a primary flow, this method 
of flow augmentation is not to be used here because of the resulting 
unsteady motion of the jet [7]. Rather, the injected flow is augmented 
as described in Section IV. 

Reference [8] states that the jet profile may be determined by con- 
Sidering the equation of conservation of charges. The radius of the gas 
jet (r) may be described as a function of the initial radius ee 
initial charge density (0 44): the applied electric field and the axial 


coordinate (z). 


he at aw leita —; 2 (5) 
+e 
Vo KE app] 


This is based on the assumption that the charge density is a function only 
of axial distance from the point of injection and that the particle 


velocity due to coupling with the neutral gas and applied field is much 


He 





greater than that due to the space chrage field, xen 

Vv tkE._ o> kE (6) 
Neglecting space charge effects disallows any radial velocity. Yet 
there must be some component of velocity. For flow along a streamline, 


Crocco's theorem for a steady state, which relates enthalpy (ho) and 


entropy (S) gradients to vorticity, is given by: 


Vx VxVS vh, - TVS 
where Vx Vv. =w is the vorticity. For the streamline defining the 
boundary of the gas jet, such gradients will be very large, giving rise 
to vorticity. The radial velocity component introduced will give a 
larger jet profile than that predicted by Eq. (5), leading to a reduced 
charge dens i ty at any given axial distance. 

‘In general, Equations (2) and (4) are coupled through the p, term, 


and must be solved simultaneously. Expanding the charge conservation 


equation: 


For particles of low mobility, and fields below breakdown, v_ will 
predominate, and Eq. (6) will hold. If the jet is assumed to be flowing 


normal to two plane parallel electrodes: 


pe Oe a ee 
ar “py ~ ar Yo * K 36 Happ 


0 


Then an assumption of an initial homogeneous distribution of charge at 


the jet entry plane will remain true at any position downstream. Thus, 


or Pe 0 
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and the charge conservation equation becomes: 


as 2 Pp) — i a Zs 
V- (o,¥)) = Pe 97 “p + 7 aane 0 


This allows the cloud shape and hence the charge density to be determined 
solely by the fluid dynamics of the dielectric jet, and uncouples the 
charge conservation equation from Poieaones equation. Another approach 
to this problem is to determine the jet profile experimentally. An 
approximate analytical description is obtained by fitting a parabola to 
this profile. This approach is convenient in this work and the 
approximation iPro an integral part of the computer program. 

We can now discuss the assumption of Eq. (6) that the charged par- 
ticles are affected only by the gas jet and the applied field. When 
charged particles are injected into the conversion region, they will face 
the field produced by the charged particles already in the conversion 
section. Eventually, these charged particles will face fields which may 
repel all charges but those with the highest degree of viscous coupling 
to the gas flow. In this manner, space charge effects provide a velocity 
"filter" through which the gas jet must move [9]. 

The expression for current density (J) in the dielectric jet is 
given by: 

T= 06 V, 
This expression would require a greater charge density as space charge 
effects slow down the charged particles in order to maintain a given 
current density. But as mentioned previously on Page 11, breakdown of the 
dielectric will occur if the space charge electric field strength is | 
greater than the dielectric breakdown strength. This tradeoff limits 


the performance of the EGD generator. 
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III. PROGRAMMING CONSIDERATIONS 


The computer program was written for the numerical solution of Eq. 
(4), since closed form solutions to Poisson's equation are difficult to 
obtain, except in very simple cases. The program is capable of solving 
Poisson's equation for either an axisymmetric or two-dimensional geometry. 

To facilitate the solution, Poisson's equation was normalized as 
described in Appendix A. In the axisymmetric form, the equation to be 


solved becomes: 


’ 
Perk | ORRe gee a (7) 


where A is the normalized potential and R is the normalized radius. The 
term on the right side is a result of the normalization described in 
Appendix A. In the two-dimensional equation, there is a change of 


variables and the radius term drops out to give the form: 
me Ye 


By covering the geometry to be studied with a uniform Square mesh, 
the problem may be defined in terms of a finite number of discrete points, 
at each of which the solution to Eq. (7) is desired. The boundaries of 
the mesh are determined either by the presence of electrodes at a con- 
Stant potential, or by a zero normal derivative of the potential $9. 
The zero normal derivative is found on an axis of symmetry, on a boundary 
with no normal current flow, such as an insulator, or on a boundary 
between two electrodes in the absence of space charge and end effects. 
End effects result in a distortion of the distributions about the end of 
a pseudo-infinite plate because of the radius of curvature of the end. 


Figure (2) depicts a sample geometry with appropriate boundary conditions. 
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The derivation of the finite difference equations for the two forms 
of Poisson's equation is given in Appendix 8B. After solving for the 
potential at a particular point, the form of the equation is: 


wile HL ; 
5 i Ty ete gt ee A 51) 


— 


where the subscripts represent mesh row and column, respectively. For 
points governed by a boundary condition of a zero normal derivative, as 


on an axis of symmetry, the general equation becomes: 


A. J = 1.333 Ay 541 a souS AG 542 


The derivation of this equation is also given in Appendix B. 

The solution method chosen is an iterative process known as "success- 
ive over-relaxation" (SOR), or the "extrapolated Liebmann method [10]." 
Over-relaxation introduces a factor (w), by which the solution at the 
point in question, as a result of the previous iteration, is taken into 
account. By assigning the proper value to this factor, convergence rate 


of the solution may be maximized; however, w is limited by the following: 
1 sms 2 


The method of selecting a value for w may be found in textbooks covering 
iterative solution methods such as Reference [10]. For this thesis w is 


a constant equal to 1.5. The form of the resulting equation then 


becomes : 
n+] _ an oy, ne n n+] n n+] 
my ti 4 G wade i Riis Oc 


nies an _ ant] n | 
+e (AP 5 AG jl) 4 s| 


7 





This shows that the solution at a point, on the (n+1)st iteration is 
dependent on the solutions at surrounding points obtained on either the 
(n)th or (ntl)st iteration. Convergence is determined by comparing a 
specified value, DELU, to the maximum potential change per iteration. 
The value of DELU is a data item required by the computer program and 
remains constant until changed by the user. 

This simplification of the governing equations introduced in 
Section II allows the space charge density to be described as a function 
of the initial space charge, of jet radius and of the initial jet radius. 
The expression, as shown by Minardi, is: 

Teg (8) 
This leads to a charge density varying with the longitudinal coordinate 
only and a constant distribution in the plane normal to the dielectric 
flow. The use of a Known nozzle radius, i.e., intial jet radius, of a 
gas flow speed, and of a corona current (I,), allows determination of 
the injected space charge density. 

Also to be considered when modeling the charge distribution is the 
collection process. Charged particles cannot continue to be collected 
indefinitely. Eventually, at some point downstream from the initial 
point of collection, all but a negligible part of the charged particles 
that are to be collected, will have been collected. Further charge 
contributions to the collector current will have a negligible effect on 
the potential distribution. It is. necessary to determine where this 
takes place, as the boundary conditions on the problem include a zero 
charge density. Hence, the space charge contained in the HENeCUnTeee 
is not to be carried up to the boundary of the region of interest. 


Using various collection lengths in the computer program, it was found 
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that allowing space charge to exist over 80% of the collector length 
gave results which changed little, less than 3%, as this collection 
length was increased. 

In general, the jet profile does not coincide with a mesh point at 
any given axial coordinate. However, the program is not capable of 
handling a variable mesh size. In order to differentiate clearly 
between the regions of space charge and no space charge, the program 
approximates the gas jet with a series of steps in order to meet the 


constant mesh interval limitation. This is explained with the help of 


the figure below: 


= 59S bl 


rtar 
eee 
1 58) 3 »J 


If ar > , the gas jet is considered to be passing through the point 


— 


(i,j+1), and if ar < om, through the point (i,j). The charge density is 
computed according to Eq. (8), giving a value of charge density which is 
then assigned to each mesh point out to, and including, the point on the 
approximated jet boundary. This introduces some error into the calcula- 
tion process. The approximation scheme limits the displacement of the: 
dielectric jet boundary to a maximum of HL/2. Then the maximum error 


(o 


ae) is given by: 


oo HL/2R 
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This would imply that the maximum error would occur at the minimum jet 
radius, the injector nozzle. By selecting the appropriate mesh interval, 
the radius at the nozzle can be made to coincide with a mesh point. The 
minimum radius at which the maximum displacement will occur is then 
downstream from the nozzle. Cie ae 5% is typical of the various mesh 
sizes used. 

The validity and accuracy of the program was determined by compari- 
son with the known solution for several problems. A geometry consisting 
of finite parallel plane electrodes, approximating infinite parallel 
plane electrodes, with a gap spacing of 2 cm and a potential difference 
of 15 kV was used to determine if the program would predict the exist- 
ence of electrical breakdown across the gap as specified in [11]. 

Because of the method of normalization used, as described in Appendix A, 
the result of entering this configuration into the program is a potential 
distribution which varies linearly from zero to one, and a field strength 
of one throughout the gap, indicating that breakdown of the air gap is 
predicted. To attain this solution, it was found that electrode end 
effects had to be minimized, or a very distorted distribution would 
result. A pseudo-infinite plane electrode minimum length to air gap 
ratio of one was required to minimize the distortion. Figure (3) depicts 
the correct potential distribution calculated by the program once, the 
parallel plane electrodes were properly modeled. 

Results of the program were also compared with two results previously 
published by Minardi [2,8], concerning space charge flow between parallel 
infinite plane wiser rode values of the input data used were the same. 

A comparison of the graphical output appeared to be equivalent for the 


two geometries compared. 


20 





IV. DESCRIPTION OF APPARATUS 


The need to input certain parameters to the program dictated the type 
and extent of the experimental work to be accomplished. This work also 
had to be relevant to the EGD facility that was under operation at NPS 
[4,5]. The program requires both the electrode geometry of interest and 
the dielectric gas jet profile, as inputs, so that charge density distri- 
bution may be calculated. The gas jet profile is also the means of 
modeling the change in a variable such as the dielectric flow speed. The 
effect of an increase in v_ was of particular interest because a larger 
v,, would increase the dominance of the first term on the right of Eq. (1). 
Also of interest was the effect of the guard ring electrode on the space 
charge flow. These two experiments, increasing Min and the use of the 
guard ring, were also required for correlation of the trends in experimental 
results with those trends predicted by the computer program. 

The EGD Bahinallkca be indeed consists of a charged aerosol injector 
and a collector probe. See Figure (4) for dumansinee The vapor for the 
aerosol is generated by a modified Hotshot model MB-31 electric steam 
boiler, pressure fed with distilled water. To obtain best results, it 
was found necessary to ensure that the entire steam generating system 
was free of any impurities. The corona needle is grounded through a 
Simpson microammeter, while the attractor ring is powered by a Sorenson 
high voltage supply, with a range of zero to ten kV. The power supply 
output is monitored by a Sensitive Research high impedance voltmeter 
connected across the corona attractor ring gap. The collector needle 
is grounded through a Simpson microammeter in series with a diode to 


prevent reverse current flow [7]. 
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In order to determine the jet profile, Schlieren photographs of the 
steam jet were taken. The setup consists of mirrors, a collimated light 
source, a knife-edge and a camera. Photographs of the jet were taken 
with stagnation pressures of 8, 10 and 14 psig. A horizontal knife-edge 
was used to more clearly define the vertical gradient across the jet 
boundary. The Schlieren photographs of the 8 and 14 psig jet are shown 
in Figure (5). The photographs show a high pressure jet (14 psig) 
approximately 100% larger than the given by Eq. (5). The resulting jet 
profiles are shown in Figure (6). 

In order to determine the effects of manipulating the electric field 
distribution, a field entirely separate from that of the corona was 
applied by means of a ‘guard electrode’, in the form of a ring of stain- 
less steel. Internal diameters used ranged from 3/8 of an inch to 3/16 
of an inch. The ring was mounted on a plexiglass insulator, supported 
by a stainless steel rod. The rod was mounted on a vernier allowing the 
position of the ring to be varied along the longitudinal axis of the 
collector needle. The voltage applied to the guard ring was supplied 
by a Spellman high voltage power supply. The output voltage of the power 
supply was determined by calibration of the power supply voltmeter with 
the Sensitive Research voltmeter. This was necessary because the 
Spellman voltmeter does not indicate line voltage. The guard ring 


Schematic is shown in Figure (7). 
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V. EXPERIMENTAL PROCEDURE 


Initial experiments with the guard ring were carried out at eight 
psig. The temperatures used corresponded to vapor states at either 
side of saturated vapor, similar to previous EGD work [4,5,7]. At this 
low pressure, a severe deterioration of the ‘convection current’ ratio, 
(collector current/corona current), was encountered as ring voltage was 
increased. This appeared to be the result of too low a dielectric 
velocity and hence a high slip parameter (i.e., a high ratio of the drift 
velocity to the dielectric velocity). The charged particles were diverted 
from the collector, being attracted to the guard ring. For this reason 
and for the desire Eee the generator performance without the ring, 
it was decided to increase the flow rate of the dielectric by increasing 
the reservoir sieSiike Stagnation conditions were changed by the 
introduction of a separate ‘dry' gas into the steam boiler. 

The baseline characteristics of the EGD channel were determined by 
Supplying steam to the corona discharge at eight psig and 240-245°F. 
Corona voltage was varied over a range of zero to three kV, and values 
of corona current and collector current were recorded. The collector 
probe was positioned 5 mm downstream from the nozzle exit plane. Dry 
gases used were both No and A. The choice of gas appeared to make no 
difference in the results. Introducing the dry gases was done by 
regulating the flow of gas from a standard high pressure gas bottle. 

The gas was introduced into the boiler through a-check valve operating 
at 10 psig. With the boiler operating at 8 psig, the high pressure gas 


was used to raise the operating pressure of the boiler to approximately 
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13.5 psig. After allowing the rate of steam production to stabilize, 
resultant stagnation pressure was 14 petite 

To determine the effects of the ‘guard electrode’, a baseline 
performance with no guard ring was run at 14 psig. Once this was 
established, the trailing edge of the guard ring was placed in one of 
three positions, either 2 mm upstream, coincident with, or 2 mm down- 
stream of the collector needle tip. At each guard ring position, the 
applied potential was varied over a range of 0 - 3.5 kV, and effects on 
collector current were recorded. Rings of various sizes were used. 
When rings smaller than 3/8 of an inch diameter were used, breakdown 
occurred, apparently caused by stray droplets providing a direct path 
for current transmission. Condensation of the steam inside the ring 


was also observed, and was a possible factor contributing to breakdown. 
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VI. RESULTS AND DISCUSSION 


The introduction of a high pressure gas into the system resulted in 
collected currents of two to three times greater than the low pressure 
steam. Figure (8) shows the baseline (8 psig) corona current-voltage 
curve and the corona current-voltage curve for the high pressure steam 
jet. The two curves are representative of the range of corona operation. 
Figure (8) also shows the variation of the value of the ‘convection 
current! ratio. 


In the expression: 


the term Ve may be given by the following expression from Shapiro [12]. 





p 2 
where: om 
lv. 
M= 
a 


fae] 


For the pressures considered, 8 psig and 14 psig, the respective Mach 
numbers are 0.832 and 1.02. These correspond to flow speeds, at the 
injector nozzle exit plane, of 402 m/s at ee = 22.7 psia, and 483 m/s 
at a = 28.7 psia. From the Schlieren photographs the profile of the 
jet was determined. Assuming continuity, this leads to a prediction of 
the variation of flow speed in the flow direction. 

The increased effectiveness of the high pressure jet may be shown 
by recalling the slip parameter, 6, defined to be the ratio of the par- 
ticle drift velocity, kE, to the steam jet flow speed. In terms of the 


non-dimensional electric field output of the computer program, E 
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the slip parameter is: 7 
| KE, [E| 


oe = eae 
Iv.| 


oO 





For a jet with particles well coupled to the flow, 6 will be less than 
one. To determine 6, the value of E may be solved for by the computer 
program. From Eq. (8), the value of Pag» the initial charge density, may 
be Pe caanined: For the low and high pressure jets, this results in 
0.0102 and 0.0084 C/m> respectively. Using this density as the program 
input parameter RHOZRO, an le | field distribution may be obtained. This 
was done for both low and high pressure cases. OQther significant para- 
meters used were attractor voltage of 2300 volts, and breakdown field 
strength of 3 x 10° V/m. Figures (9) and (10) are plots of the equi- 
potential surfaces resulting from the computer solution of the low and 
high pressure jets. Using these results, the variation of 6 along the 
centerline of the EGD channel was obtained, and is plotted in Figure (11). 
The values of 6 are based on an assumed mobility of 4 x 10°? m@/V-s, 

a value representative of the range of mobilities of the EGD charged 
particle [7]. 

Figure (11) shows that everywhere along the centerline of the EGD 
channel, with the exception of a single point, the high pressure jet has 
a lower slip than the low pressure jet. The single point that differs 
may well be the result of the inaccuracies of the numerical solution. 
The use of a higher value of mobility would tend to give values of 6 
greater than one. However, if the conjectured distribution of [5] is 
approximately correct, particles of higher mobility are a small part of 


the total number of charged particles. 
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From Eq. (3), the total E field has two contributions. One, Fone 
is due entirely to the geometry and voltages applied to various parts 

of the EGD geometry. This can be changed only by a physical change in 
the geometry or applied voltages. The potential and E field distribution 
resulting from the geometry alone may be solved for by Laplace's equation. 
The Laplacian solution represents the potential or E field distribution 
unaffected by the presence of space charge. Figure (10) is the solution 
to Poisson's equation, and it shows the added distortion of the dis- 
tribution due to the presence of space charge, when compared to the 
Laplacian solution, Figure (12). The added distortion appears in Eq. 


(3) as the second contribution, E.., the space charge field. Figure (13) 


SC? 
depicts the \E| field distribution corresponding to the potential 
distribution of Figure (10). Since the space charge is of the same 
polarity as the charged particles, their movement is opposed. Elimina- 
tion or reduction of : would enhance charged particle velocity, 
making the collection easier. Application of a potential to the guard 
ring is designed to impose an additional external field on the flow. 
However, if the field HVE LCOMpensates for the space charge field, and 
the ring becomes a dominant sink for the charged particles, the purpose 
of the ring is defeated. 

Figure (14) shows the increase in collector current as a result of 
using the guard ring. ie results are shown over a range of voltages 
at the three positions used. It appears that the position achieving 
the greatest increase is that of the trailing edge plane coincident with 
the tip of the collector needle. From Figure (10), this is also the 
area of greatest space charge distortion, implying that the effects of 


the ring are concentrated in a small area. 
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Figures (15), (16) and (17) are graphical computer outputs of 
potential distribution for a grounded ring, and potentials of | and 2 kV. 
The experimental results also show the proper trend in collector current 
as ring voltage is increased. That is, the current slowly increases to 
a broad maximum, and then drops off with further voltage increases, the 
maximum occuring at about | kV. 

The computer predicted results of using the guard ring is shown in 
Figure (18) as a plot of centerline Laplacian E minus the Poisson E 
field. As this difference approaches zero, the particular guard ring 
distribution approaches the optimum Laplacian distribution. The computer 
results indicate that the E field is definitely modified by the presence 
of the guard ring. The predicted optimum ring voltage, corresponding to 
maximum collector current, appears to lie between ground and | kV. The 
point common to all curves lies at the midpoint of the ring, and is prob- 
ably due to the uni formi ty of the guard ring electric field at this 


point. 
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VII. CONCLUSIONS 


The computer program appears to present an accurate picture of the 
potential and electric field distributions within an EGD device. This 
is borne out by the computer prediction of experimentally observed trends 
in EGD channel performance. The ‘convection current’ ratio, I/Iy; was 
found to be affected by the steam reservoir conditions. In addition, 
computer and experimental results agree and show that increased EGD 
performance can be attained by suitable manipulation of the electric 
field distribution, and hence of the space charge flow. 

The agreement found between analytical and experimental results 
Suggests that the use of the computer program as a design tool, par- 
ticularly when modeling the space charge flow in a high pressure jet, 
is feasible. 

To be completely general in approach, some further account should 
be taken of particle mobility. Since the drift velocity appears to be 
negligible at high pressures, this is especially true when working with 
particles of high mobility, and/or dielectric jets at low speeds. In 
this case, the slip parameter could go well above 1, and the present 
assumption of zero slip will give completely invalid answers. Some 
preliminary work on determining the mobility range of the charged 


particles used in these experiments is shown in Appendix D. 
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VIII. RECOMMENDATIONS 


In addition to the measurement of mobility, the difference equations 
should be modified to allow for a variable mesh size. Equations of this 
type may be found in a textbook such as Reference [10]. This will allow 
a greater generality in selecting geometries to be studied, and will also 
do away with a great deal of the error involved in discrete jet profile 
approximation. It vail almost certainly mean greater computer storage 
requirements. 

As a design tool, the program can be used to determine the optimum 
guard ring geometry and voltage. From the curves of Figure (18), the 
solutions for the various voltages diverge from the optimum in the 
vicinity of the collector needle. A guard ring geometry capable of 
affecting this divergence, yet maintaining the distribution of the 
remainder of the 1 kV curve, would appear to be an improvement. 

It should also be possible to use the program as a design tool to 
optimize both injector nozzle and collector geometries. A nozzle 
geometry allowing for the use of smaller guard ring sizes may allow 
for more effective guard ring performance, yet forestall breakdown. 

An optimum collector geometry may exist which will allow for greater 
Charge collection efficiency resulting simply from the geometry, rather 
than by flow manipulation to achieve the higher collector currents. 

Other applications outside the EGD field may also exist. For 
instance, the program may be used to study the region of space charge 
known as an electrode 'sheath' which Ravidee a neutral plasma and an 
electrode. Also, Reference [4] covers a study of EGD control of separated 
flow. This study may be extended through use of the program to determine 


optimum electrode geometries for control. 
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The vector form for Poisson's equation is: — 


After expanding in cylindrical coordinates: 


By defining: 


where: 
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The normalized equation becomes: 
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By combining and clearing the left side: 





Similarly for the two-dimensional form of Poisson's equation. 


defining: 
as Saye 
X h : i . 
The two dimensional form becomes: ms 


The electric field may be normalized by defining: 


2 a 
E = es 
ED 
where: E = the normalized value of E 
EY = the dielectric breakdown field strength 
from: E = -Vo 
HE. _ shige oe. 
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APPENDIX B: FINITE DIFFERENCE EQUATION DERIVATION 
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Knowing the value at any point in the mesh, the value at another 





point, nearby, may be found by means of a Taylor series. By writing the 


Bemres in the J direction, expressions for A. and A. El may be 


ait 
found. These are: 
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Aa 34] = AS + HLA; + oT Aaj +m Ay ae (1) 
A = A = AEA. . + HL ai - HLT Ag .r (2) 
1,Jj-] eg] Melee - lis) “Os ~ Wee” wae 
Solving for A. rE neglecting higher order terms where primes indicate 
the derivative with respect to the direction of interest: 
Al w= A, Jt] ‘- J-] (3) 
ie): 2HL . 


To determine the finite difference equation for the second partial 
derivative, add (1) and (2). This results in: 


~ 2 it 
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Solving for As ye again neglecting higher order terms: 


isi, ~ sell (4) 


By writing the series in the I direction, solving for values of A. 4) j 


and A. 4 5? and adding the resulting equations, the second derivative 


in the I direction may be determined. The resulting equation is: 


Cea ee eee 
Be ede le (5) 
sd HL 


By combining (3), (4), and (5), and solving for A. , a finite difference 


oJ 
form of Poisson's equation in axisymmetric form is obtained. The 


axisymmetric form is: 


BR. . = Ui co em 
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For points along a boundary where the normal derivative is equal to 

zero, the equations for the value on the boundary may be derived in the 

Same manner and cast in the form of second order forward or backward 

difference equations. Assuming (1,j-1) is on the boundary of interest, 


the series equations become: 


2 
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Multiplying (8) by 4 and subtracting from (7): 
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Neglecting higher order terms: 


- 4A 3 A 
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Setting the normal derivative equal to zero, and solving for A. j-1° 


A = 1.333 A. 57 som (10) 


1 sJ-| 1 sJtl 
The successive over-relaxation method (SOR) uses previously calculated 
values to speed up convergence. For a point inside the mesh, the final 


difference equation becomes: 


ae lov cowall Ww) 2 n n+ n n+] 
Pye ATS TH [MPC + aS AT tA gent ATI 
(11) 
HL n n+] n 
* an (i,j 7 ijt) - ae 
For a point on a boundary of zero normal derivative, the equation 
becomes: 
n+] _ an n n n 
A. “el ae A. ica + oft38 A. J = 2333 A. jt = A. a (12) 
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APPENDIX C: PROGRAM DESCRIPTION 


The computer program is designed to solve a finite difference form 
of Poisson's equation, either in axisymmetric or two-dimensional form, 
at each node defined by a uniform square mesh covering the geometry to 
be studied. In order to do this, the program must be able to do the 
following: 

1. Determine which finite difference equation is to be solved. 
This is : function of the geometry and the node position within 
the mesh. 

Determine and store charge distribution at each node. 
Solve the difference equation and store the answer. 


Determine when a solution to the problem has been attained. 


no > WW P> 


Provide the results in an easily interpreted form. 

The program is made up of two major parts:. preparation and calcula- 
tion, and output of results. In general at each point in the EGD channel 
defined by the uniform mesh, there exists a: charge density and a potential. 
The charge density may be zero. This can be the case inside of the jet, 
as well as outside. The voltage may either be the result of an electrode 
presence or a potential surface between electrodes. If the voltage is 
due to an electrode presence, this value will remain constant, otherwise 
it will vary as a function of the charge density. The program must be 
able to make this differentiation. In all, the program must know the 
potential at a point, the charge density, and whether the point aerine: 

a portion of an electrode. In order to simplify the program coding, the 
program stores all values of potential in one array, "A". The values of 


charge density are stored in a separate array, "G", of the same size. 
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This array is also used to store all information required to determine 
whether the potential is that of an electrode, the potential description 
value (PDV), or is to be changed as a function of the charge density. 

This can be done since a point of constant potential will be an electrode, 
and a charge density cannot exist at the same point. 

The program determines convergence of the solution method by compar- 
ing answers from the (nt1)st iteration with those obtained on the (n)th 
iteration. This will reveal the maximum change in potential as a result 
of the iteration. If the maximum is less than some designated value, the 
iterative procedure is stopped. This method requires storage for the 
answers generated by the (n)th iteration. These are stored in array "B". 
The value of the potential resulting from the previous solution is required 
for the SOR method. This value is supplied from array "B". Hence the 
variables required for the solution of the finite difference equation 
are taken from three separate sources. A superposition of arrays "A" 
and "B" describes the complete situation at any node within the EGD 
channel. 

The program has been written for use on an IBM 360/67 conliter 
System. The output is in both printed and plotted form. The plotting 
was accomplished on the CALCOMP Model 765. The program consists of 


the following subroutines: 


SUBROUTINE NAME PURPOSE 

MAIN (Main Program) General coordination of preparation 
and calculation phase and output 
phase. 

GSPRD Calculates cloud radius. Determines 


charge density at each point, and 
Stores in "G" array. 
GENFIL Reads geometry input. : 
CALC] Calculates potential. Determines E 
field components and modulus. 
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WRITER Prints values of potential, E field 
and E field components, determined 


in CALC]. 
CLVL Determines lines of constant potential 
and constant E field to be plotted. 
FLOP Orients arrays to coincide with output 
labeling convention. 
OUTPUT General coordination of output phase. 
CONTUR Coordination of labeling and plotting 


of graph, and preparation for plotting 
results. Modified NPGS Computer 
Library subroutine, which allows the 
Superposition of three sets of contours 
on one plotted output. 

All aay except the electrode geometry is read in under NAMELIST 
control. This requires a special format for all cards, but provides for 
a more flexible means of handling the data in that there are no fixed 
fields in which the variables must appear, nor must the arrangement 
coincide with that of the NAMELIST statement. 

The first data card required by MAIN is the card denoting the appear- 
ance of the data group. It has the format, starting in column 2; &PARAM. 
PARAM is the NAMELIST name. The first column of all data cards is ignored. 
Each data item may be separated from the following by a comma. However, 
no comma must appear between the NAMELIST name and first data item. 

More than one data item may appear on a card. Each data item is in one 
of the two forms below. 

1. “Variable name" = FORTRAN constant 

2. "Array name" = a set of FORTRAN constants, separated by constants 

Lis). 
The variables read in are listed below. The first five are explained 


with the help of the representative geometry in Figure ( 2). 
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NAME TYPE DESCRIPTION 


NATRAD 1*4 Mesh column number of attractor inner 
radius. This is the point defined as 
(6,3). Then the value assigned to 
NATRAD is 3. 

NBLOCK 1*4 The number of discrete points, lines, 
or blocks of points required to 
completely describe the charge dis- 
tribution or geometry of the problem. 
For the figure shown, the attractor 
ring may be defined by a block of 
points, a line and a single point. 
The corona needle, and collector 
needle may each be defined by two 
lines of points. In addition, if it 
is assumed that the collector needle 
collects all the charge, a block of 
points can be used to define the area 
within the dielectric jet where in no 
charge exists. Thus, for the geometry 
and charge distribution of the problem 
Shown, NBLOCK = 8. 


NDLPNT [*4 Mesh row position of corona needle 
tip. For the geometry shown, NDLPNT = 4. 
NDLRAD [*4 | Mesh column number of corona needle 
| centerline. For geometry shown, 
NDLRAD = 1. | 
NOZPOS  «T*4 Mesh row number of attractor exit 
' plane. For the geometry shown, 
NOZPOS = 6. 
BRKDWN R*4 Breakdown field strength of dielectric 
gas. (volts/meter) 
CONVER R*4 Characteristic channel length of EGD 
geometry. (centimeters ) 
DELU | R*4 Value compared with maximum difference 


between iterations. If the maximum 
difference is less than desired value 
of DELU, convergence of the iterative 
Scheme is signaled. 

GAPPOT R*4 Maximum fixed potential covered by the 
mesh. In general, the attractor will 
be at the maximum potential. 


TH [*4 Desired height of output plot measured 
along longitudinal coordinate (IH < 15) 

IW [*4 Desired width of output plot measured 
‘along radial coordinate (IW <9) . 

NCOLS [*4 Number of mesh columns (NCOLS < 61). 

NDIMEN = 0 [*4 Indicates problem is a two-dimensional 
geometry in constant area channel. 

= | [*4 Indicates problem is axisymmetric. 
NGLOSS = 0 [*4 No charge density decrease accounted 


for in the conversion section other 
than by charge cloud spreading. 
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NGLOSS = 1 I*4 Some loss mechanism (recombination, 
| etc.) and cloud spreading will be 
accounted for in determining charge 
cloud densities. 
NL [*4 Number of potential contours desired 
in plotted output (NL < 20). 
NLI [*4 Number ‘of E field contours desired in 
plotted output (NLI s 20). 
NOLINE = -1 [*4 Only potential and E field contours 
will be plotted. 
= 0 1*4 Only physical geometry (fixed potential) 
points will be plotted. 
= | [*4 Physical geometry will be superimposed 
on contour plot. 
NROWS I*4 Number of mesh rows (NROWS < 81) 
RHOZRO R*4 Charge density at exit plane of 


attractor (coulombs/cubic meter) 

The last card in the NAMELIST signals the end of the list. The format 
is &END. Again the symbol ampersand must appear in column 2 of the data 
card. 

The jet radius and charge density distribution are calculated by 
Subroutine GSPRD according to the charge profile given by statement 
numbered 35. In the case of a two-dimensional geometry, only the charge 
distribution is calculated, as the jet is assumed to fill the entire 
two-dimensional channel. The charge density will remain constant through- 
out the flow in the channel, unless some charge loss mechanism is 
introduced, signaled by the value of NGLOSS. The jet profile and charge 
distribution for a geometry with an off center corona needle is also 
determined by GSPRD. 

Subroutine GENFIL is used to input the specific values of potential 
and charge density or PDV. This is accomplished with the help of the 
variable NBLOCK, which is used as a DO-LOOP parameter. The loop causes 
the program to read the exact number of data cards defined by NBLOCK. 
Since each block defined by NBLOCK contains points of common value, the 
entire block can be represented by a single data card. The format for 
these data cards is 413,2F12.4 

40 





The variables read in by GENFIL are: 
NAME TYPE DESCRIPTION 


NII [*4 Mesh row nunber of beginning of line 
or block. 

NI2 I*4 Mesh row number of end of line or 
block. 

NJl [*4 Mesh column number of beginning of 
line or block. 

NJ2 1*4 Mesh column number of end of line 
or block. 

AV R*4 Normalized potential common to all 


points encompassed by NII, NI2, 
NJI and NJd2. 

GV R*4 EITHER: Normalized charge density 
common to all mesh points encompassed 
by NIT, NI2, NJI and NJ2. 

OR: Potential description value (PDV) 

if the defined points are at a fixed 

potential. 
An arbitrary value of 99. is the PDV used to describe a point of fixed 
potential. No arithmetic operations will be carried out at this point. 
Note also that NI] < NI2, and NJl s NJ2, since these variables are the 
beginning and end of a DO-LOOP. In the case of a point, NII = NI2 and 
NJ 1 = NJ2. In the case of a line of points, either NI] = NI2 or NJl = 
NJ2, depending on the orientation of the line. The values assigned to 


the above variables to describe the representative geometry are shown 


below: 
NI] NI2 NJl NJ2 AV GC - Corresponding to 
] 6 5 15 ea 99, NBLOCK 1 
5 6 4 4 1.0 99. NBLOCK 2 
6 6 3 3 1.0 99. NBLOCK 3 
] 4 ] ] 0.0 99. NBLOCK 4 
] 3 2 2 0.0 99. NBLOCK 5 
14 20 ] ] -0.1 99. NBLOCK 6 
15 20 Z 2 -0.] 99. NBLOCK 7 
17 20 2 6 O25 0.0  NBLOCK 8 


At the same time that the values of AV and GV are being loaded into the 
appropriate arrays, the coordinates of those points with a PDV of 99. 


are being converted into plotter pen position coordinates to aid in 
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plotting the output requested by NOLINE. The pen position coordinates 
are stored in the appropriate vectors, XLIN and YLIN for later use during 
the output phase. 

Subroutine CALCI solves the appropriate form of the finite difference 
equations, depending on the geometry variable NDIMEN, and location of the 
point in question. In addition, once the potential distribution has been 
calculated, CALC] determines the E field, and the longitudinal and radial 
components of the E field. The value of the E field and its components 
are normalized with respect to the breakdown field strength of the 
dielectric, BRKDWN, while the values of potential are normalized with 
respect to the maximum fixed potential, GAPPOT. The derivation of the 


normalization of E is shown in Appendix A, which gives: 


ge 
E = - hE VA 


The characteristic EGD geometry length 'h' is read into the program by 
the value of the variable CONVER. Defining: 


| 
SALPH = hE, 


JE] = $ALPH |-vAl 


The value of the constant $ALPH is calculated by MAIN, and is passed to 
CALC] as a calling argument. By normalizing E in this manner, a quick 
glance at the results will show whether or not the electric fields pro- 
duced will exceed the dielectric breakdown strength. The results of 
CALC] are printed by subroutine WRITER. 

For the graphical output of the results, subroutine OUTPUT coordina- 
tes the various plotting subroutines required. The remainder of the 


required data cards are read in by subroutines OUTPUT under NAMELIST 
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control. The first data card is &TABL, starting in column 2. The 
variables in NAMELIST TABL are read into two vectors, TABLE and LTG, 
in the format: 

"Array name" = a set of FORTRAN constants, separated by a comma. 
The format for the elements of TABLE is TABLE( ) = '-------- ', where (-) 
indicates a card column available for a data item. This format prepares 
the data to be used in labeling the graphical output. The variables 


read in are: 


NAME TYPE DESCRIPTION 
TABLE (1) R*8 Nozzle radius in centimeters 
TABLE (2) R*8 Initial charge density in C/m 
WeApLe (3) R*8 Maximum fixed potential (GAPPOT) 
TABLE (4) R*8 Dielectric premittivity 
TABLE (5) _ R¥8 Breakdown field strength 

of dielectric 
TABLE (6) R*8 Characteristic length modeled 


The elements of the vector LTG can be used to construct a grid for the 
plotted output, label it and the plotted contours. The vector LTG is 
typed LOGICAL * 1. The variables read into the vector LTG are: 


LTG (1) ea ~TRUE. All exterior contour segments 
and those interior contour 
setments which represent a 
local maximum will be labeled. 

~FALSE. Omits labeling option. 

LTG (2) L*] ~TRUE. Request tic marks and corres- 
ponding scale values one inch 
apart on the exterior frame 
of the contour graph. 

SEALSEs s Omits: cb sODeron. 

G(s) L*] .~TRUE. Request a one inch by one 
inch straight line grid to 
Superimposed on the contour 
graph. 

FALSE. Omits the grid option. 


Again, the last card signals the end of the NAMELIST. The format is: 
&END 


with the symbol ampersand in column 2 of the data card. 
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APPENDIX D: MOBILITY MEASUREMENTS 


The apparatus used to investigate the mobility of the charged particles 
is shown in Figure (19). This consists of a fine Monel steel grid con- 
centric with a stainless steel cylinder. The EGD process under investiga- 
tion involves the transport of negatively charged particles. Hence if 
a voltage is applied to the grid and varies from zero to some value below 
ground, the grid will repel or pass negatively charged particles. If at 
the same time, the outer cylinder is maintained at ground, those particles 
passed by the grid will be accelerated across the gap, giving a noticeable 
variation in current flow due to the charge arriving at the outer cylinder. 
From the definition of drift velocity: 

v= kE 
and assuming that the E field is uniform: 


ieee ae 
v =k ar ais 


where | is the spacing between the grid and outer cylinder. But since: 


Vet 

The mobility may be determined by: 
ee 
tV 


By varying the frequency with which the grid voltage changes and ensuring 
that the waveform is a Square wave, the entire range of mobility sGeeeeeea 
by the charged particles can be estimated, assuming that mobility is | 
proportional to the electric field between the cylinders. This is 


because at the lower frequencies, all particles will be affected, while 
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at the higher frequencies, or small t, inertia effects will deter all 
but the smallest or most highly charged particles from completing their 
journey. This is similar to the method of Rutherford as outlined by 
Darrow [14]. 

The measurements leading to an approximation of the mobility distri- 
bution were made at a reservoir pressure of eight psig, over a wide range 
of temperatures, 238 - 296°F. After steam flow conditions had stabil zed, 
i1.€., proper pressure and constant temperature, the device was inserted 
into the EGD channel to such a position that it was axially symmetric 
with the steam jet, but not touching. This was done to keep condensation 
on the grid to a minimum. The grid was powered by a square wave genera- 
tor while the outer ring was maintained at ground. See Figure (18) for 
a schematic of the setup. Current from the ring was fed to an oscillo- 
scope for observation. It was found that the ring current varied con- 
Siderably as the charged particles arrived from the Monel grid. This 
variation existed over a range of fequencies. At high temperatures, 
296°F, the range of current excitation appeared to be concentrated at 
a single point, 2.1 kHz. As the temperature of the steam was decreased, 
the range increased considerably, until at 238°F, the current excitation 
existed from an upper frequency of 5 kHz to a lower frequency of.5 kHz. 
At 296°F, the mobility corresponding to 2.1 kHz is 4.15 x 107* 
me/V-s, while at 238°F, the mobility range runs from 1 x 107° m?/V-s to 
lhe x ine m@/V-s. The range of mobilities appears to arise from both the 
ability of large particles to move at low frequencies and the lower tem- 
perature producing a greater distribution in particle size. ‘However, 
further work remains to be done to determine both size distribution and 


further definition of the mobility distribution. 
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ROW INDEX (1) 


Sample EGD Geometry and Boundary Conditions. 


Figure 2. 








Figure 3. Computer Prediction of Parallel Electrode Normalized 
Potential Distribution at Breakdown. 
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Figure 4. Aerosol Injector and Collector Probe. 
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Figure 11. Computer Prediction of Centerline Distribution 
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